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ABSTRACT 

GJ 581d is a potentially habitable super-Earth in the multiple system of 
exoplanets orbiting a nearby M dwarf. We investigate this planet's long-term 
dynamics, with an emphasis on its probable final rotation states acquired via 
tidal interaction with the host. The published radial velocities for the star are 
re-analysed with a benchmark planet detection algorithm, to confirm that there 
is no evidence for the recently proposed two additional planets (f and g). Limiting 
the scope to the four originally detected planets, we assess the dynamical stability 
of the system and find bounded chaos in the orbital motion. For the planet d , the 
characteristic Lyapunov time is 38 yr. Long-term numerical integration reveals 
that the system of four planets is stable, with the eccentricity of the planet d 
changing quasi-periodically in a tight range around 0.27, and with its semimajor 
axis varying only a little. The spin-orbit interaction of GJ 581d with its host 
star is dominated by the tides exerted by the star on the planet. We model 
this interaction, assuming a terrestrial composition of the mantle. Besides the 
customarily included secular parts of the triaxiality-caused and tidal torques, we 
also include these torques' oscillating components. It turns out that, dependent 
on the mantle temperature, the planet gets trapped into the 2:1 or an even higher 
spin-orbit resonance. It is very improbable that the planet could have reached 
the 1:1 resonance. This enhances the possibility of the planet being suitable for 
sustained life. 

Subject headings: planet-star interactions — planets and satellites: dynamical 
evolution and stability — celestial mechanics — planets and satellites: detection 
- planets and satellites: individual (GJ 581) 
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Introduction 



Habitability of an emerging class of super-Earths - exoplanets with masses greater than 
that of the Earth but lower than that of Uranus - depends on a combination of physical 
parameters. Among these, crucial is the intensity of irradiation from the host star. A 
favourable rate of irradiation permits water to be available in the liquid form. The irradiation 
intensity depends on the luminosity of the star, the size of the planet's orbit and, to a lesser 
degree, on the orbital eccentricity. The chemical composition of the planet's atmosphere 
influences the average temperature and the temperature variations on the surface. For 
example, estimates show that a certain minimal amount of CO2 in the atmosphere of GJ 
58 Id is required to keep water above the freezing p oint on the surface, wh i le GJ 581c is likel y 
to have experienced a runaway greenhouse event (jvon Paris et al.l l2010t iHu fc Dins 1201 ll ) . 
Planet GJ 581d, which is the m ain target of our s t udy, m ay be located on the outer edge of 
the habitable zone, according to Ivon Braun et al.l ( 120111 ). 



We begin our study with addressing the still controversial problem of composition of 
this planetary system. In Section |2j we confirm that only the four originally detected planets 
(b through e) are real and detectable at the 0.99 confidence level by our ful ly automated, fast 
grid search algorithm. The two additional planets, f and g, proposed by IVogt et al.l (120101 ) 
are not found with any combination of the published radial velocity (RV) data. 

The orbit estimation technique employed in our paper is designed, in particular, to pro- 
duce accurate and robust results for the eccentricity of each detected planet - an important 
asset for the subsequent dynamical simulations. The dynamical stability and the presence 
of chaos in the orbital motion of the four planets is investigated in Section [3] by long-term 
integrations, assuming zero inclination and coplanar orbits. The system is found to be long- 
term stable but strongly chaotic, with the eccentricity and semimajor axis of planet d little 
varying over gigayears. This allows us to investigate, in Section HJ the spin-orbit dynamics 
of GJ 581d with its orbital parameters fixed. The planet is assumed to have a terrestrial 
rheology and a near- zero obliquity. 
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2. How many planets have been detected in the GJ 581 system? 



Although this question has attracted a lot of attention in the literature, some doubts 
seem to be lingering, and a few papers have been published on subtle features of the plan- 
ets that may well be nonexistent. The first four plane ts were discovered with the same 
ins trument, HARPS at La Silla Obs ervatory: GJ 581b faonfils et alihoosh . GJ 581c and 
d (jUdry et al.l 120071 ). and GJ 581e (IMayor et al.ll2009l ). The combined series of observa- 
tions with HARPS, published in the latter paper, spanned 1570 days, while a typical single 
measurement precision was about 1 m s _1 . 

At the first stage of our investigation, we process the HARPS data with a simple bench- 
mark planet detection algorithm. This algorithm is an iterative process of optimisation, 
based on a simple grid search in the space of free parameters. The grid search is applied 
only to fitting of the nonlinear parameters, i.e., of the orbital period, eccentricity, and peri- 
astron time. The other parameters are determined by a direct Least Squares solution. At 
the first stage of fitting, the most significant sinusoidal variation in the radial velocity dat a 
is determined by a corrected Lomb-Scargle periodogram method (jLomblll976l ; IScargldll98ll ). 
Our upgrades of the method concern mainly the way the common offset of the RV points is 
treated. Instead of the traditional subtraction of the mean RV from the data prior to the peri- 
odogram computation, we fit the entire set of model functions [1, cos(27rtj/pj), sin(27ri,/pj)] 
for each trial period Pi, where tj are the times of observations. This, seemingly trivial mod- 
ification allows us to update the amplitudes of already detected planets once a new planet 
is detected, and to keep the original RV measurements unchanged throughout the data re- 
duction cycle. The mean RV is inevitably biased when a non-integer number of waves is 
present in the data. This way, subtraction of the mean RV from the data is legitimate only 
when no periodic signals are present in the data, which defeats the purpose of periodogram 
analysis. It can be proven that the periodogram value in the classical Lomb's algorithm is 
equivalent to the absolute change of the \ 2 before and after the Least Squares fit of the 
cos- and sin-terms. However, fitting these terms to the data corrupted by the subtraction 
of the mean RV gives rise to sidelobes and biases in the periodogram. The downside of 
this modification is that a direct computation of the post-fit \ 2 by Scargle's formula is no 
longer possible. Instead, the complete Least Squares solution of the matrix equation should 
be computed for each trial period. This complication, however, becomes barely noticeable, 
taken the power of the present-day computers. 

Once the period of the most prominent sinusoidal variation gets determined, a separate 
grid search is implemented to determine the best-fitting eccentrici ty. This technique is based 



on the theory of har monic decomposition of orbital motion (e.g.. iBrower fc Clemencdll961 



Konacki et al.ll2002l ). The observed radial velocity variation due to the orbital motion of a 
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single planet can be written as 



oo 



oo 



z = - C — 1 a i( e ) cos(lM) + H — 1 <M e ) sin(/M) , (1) 



where e is the eccentricity, C and H are the Thielle-Innes constants, P is the orbital period 
determined from the periodogram, M. is the mean anomaly, while oti{e) and <pi(e) are coef- 
ficients which depend only on e and define the relative magnitudes of harmonics. Although 
the coefficients can be computed via the Bessel functions of the first kind, we find it practical 
to compute them numerically using Kepler's equation at the required resolution in e between 
and 1, and to store the resulting values of the coefficients in a table. The further processing 
is then reduced to solving the overdetermined equations 



where p\ = I ai(e) and qi = I (f>i(e) y/1 — e 2 are tabulated harmonic coefficients, while the 
integer j = 1 , 2, . . . , N serves to number the available data points. Least-squares solution 
of the problem is performed for a grid of e , and results in selection of the value minimising 
the post-fit x 2 ■ The solution vector s includes three elements: s(l) is the RV offset, while 
s(2) , s(3) are the fitting coefficients from which the Thielle-Innes constants C and H can 
be derived. In parallel, a separate grid search is carried out for the phase of orbital motion, 
i.e., for the mean anomaly at the first observation. 

Finally, we have to estimate the confidence of the orbital solution. We are using an 
adaptation of the significance F-test, which was carefully tested and verified by Monte-Carlo 
simulations on a large number of real and artificial planetary systems. In this case, the 
null hypothesis for the F-test is that the data set contains a constant offset and a random, 
uncorrelated noise. Fitting of a harmonic [cos(2irtj/p i ) , sm(2irtj/pi)] to this set by the 
least-squares method will generate a \ 2 change corresponding to a random realisation of 
the F-statistic, F» , which is distributed as Fpdf • The confidence of a detection at pi is 
simply the probability of F < Fi , which is defined by the known cumulative distribution 
function F CDF . It is customary to accept a detection if the confidence is greater than 0.99, 
i.e., if the observed reduction in % 2 can happen within the null hypothesis in less than 1 
trial out of 100. The confidence of detection is computed as 
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where xi an d x\ are the value of x 2 before and after the orbital fit, respectively; while d± 
and d>2 are the corresponding numbers of degrees of freedom. The / pow multiplier is required 
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to take into account the fact that we are testing multiple periods pi and are selecting the 
one that delivers the greatest reduction in x 2 - If the probability of a single realisation of 
F to be less than a certain limit F^ m is Fcr>p(F^m), then the probability of the largest F 
among n independent trials to be less than F hra is F C DF(Fi im ) n . Thus, the / pow multiplier in 
Equation [3] is the number of independent frequencies for a given periodogram search window. 
For irregularly spaced observation times and search intervals that often exceed the Nyquist 
boundaries, evaluation of this number is nontrivial. We employ the following estimate, 
which is analogous to the Nyquist limit: / pow = ceil(AT/2/s mean ) , with s mean being the 
mean separation between the RV data points. 
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Fig. 1. — The generali sed y 2 periodogram of the GJ 581 system prior to orbital fitting based 
on the HARPS data flMayor et al.l 120091 ). The dips corresponding to the four detectable 
planets are indicated with arrows and planet names. 



Handling of multiple planet detections is another important feature of the algorithm. 
Once a planet with period P m is detected, the corresponding terms [ ^ z pi cos(2n(tj — t 0m )/P n 
J2i Qi sin(27r(tj — t 0m )/P m ) ] get permanently added to the model. Then the search for an- 
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other planet is performed outside the already detected frequencies, on the original RV data 
(which are never altered in this algorithm). With m planets already detected, the size of the 
model is 2m+l. Also, the number of degrees of freedom d\ is reduced by 5 with each detected 
planet. The coefficients C and H and the RV offset are re-adjusted each time a new planet is 
detected, which provides for a better decoupling of planets in complex systems, especially of 
planets on commensurable or resonant orbits. Fig. [1] shows the initial x 2 periodogram of the 
HARPS data prior to any planet detection. The dips corresponding to the planets detected 
from these data are marked with arrows and planet designations. The strongest signal is 
associated with planet b, which shows up first. The planets are detected in their historical 
order: b, c, d, and e. The next strongest signal in the cleaned periodogram corresponds to a 
period of 391 d. Formally, though, it should be rejected because the confidence level is only 
0.975. 

The orbital and physical parameters estimated from the fits of the four det ected planets 



are gi ven in Table [TJ Our periods are in close agreement with the estimates by iMayor et al. 



( 120091 ). but e and u (the longitude of the periastron) are markedly different . The likely 



reason is that these parameters for the planets b and e were "fixed" to in IMayor et al. 



( 120091 ). which is not advisable for multiple systems. In particular, ignoring the detectable 
eccentricity in a least-squares adjustment leaves in the periodogram side lobes and harmonics 
associated with the already detected planets and can distort the results for other planets. 
Our estimate of e for planet d is 0.27, which is lower than the 0.38 ± 0.09 obtained in Ibid. 
Our results confirm that the planets c and e may be in the 4:1 orbital resonance, and the 
closeness of our estimated u lends more credence to this possibility. 



Our results are in general agreement with the work by iTuomil ( 120111 ) who re-analysed 
the HARPS and HIRES data together using a Bayesian approach and found evidence of 
only four planets existing. However, their estimated eccentricities are consistent with for 



Table 1. Orbital parameters of the four-planet GJ 581 system (b through e). 
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all four planets. We too applied our benchmark algorithm to the combination of HARPS 
and HIRES data, introducing separate RV offsets for the two data sets. This resulted in 
detection of only two planets, b and c, but with significantly weaker signals. The planet 
d was rejected due to insufficient confidence . Fina lly, the algorithm failed completely to 
reproduce the results obtaine d by IVogt et al.l (120101) from their data. In their most recent 
update on the GJ 581 system, IVogt et al.l (120121) used the yet unpublished, much expanded 
set of RV data from the HARPS (jForveille et al.ll201ll ). They found that a fifth planet (f) 
can be detected if the eccentricity of the 66.7-day planet (d) is set to zero. Using the same 
data set, but not rejecting any data points in it, we were able to confirm this statement 
with our detection algorithm by forcing the eccentricity of the four planets b-e to zero. A 
32.1-day planet with a projected mass of 1.99Me emerges at a confidence level above 99%. It 
is likely therefore that we are dealing with the second harmonic in the signal of an eccentric 
planet, which can be confused with a separate nearly commensurate planet. A potentially 
good way of resolving this controversy is to look for the third harmonic of the 67-day signal, 
which should be generated by an eccentric planet with that period. Indeed, we find a dip in 
the x 2 periodogram at pj ~ 22.5 days after fitting out the four planets at zero eccentricity, 
but its presence can not be taken as conclusive evidence. At e = 0.27, the amplitudes of 
the first three harmonics scale as 1:0.256:0.073, hence, the expected amplitude of the third 
harmonic i s only 2 x 0.073 ~ 0.15 m s _1 , so it should drown in the observational noise. 



Vogt et al.l ( 120121 ) make a strong argument that a system of four eccentric planet can be 



dynamically unstable. We therefore set out to establish that the original four-planet system 
as determined from the HARPS measurements is dynamically viable. 



3. Chaos and stability of the orbits 



The preceding studies of the dynamical status and evolution of the plane tary system 



GJ 58 1 have been mainly concerned with verification of its physical stability. iBeust et al. 



( 120081 ) integrated the system of three planets (b, c, and d) for 10 8 yr at different inclinations 
to the line of sight (which changes the estimated plan etary masses by a factor of sin -1 i) and 
for different orbital eccentricities. iMayor et al.l (120091 ) performed similar integrations for the 
4-planet system (b through e) with an updated period of planet d. In all these simulations, 
the orbits were assumed to be coplanar. The systems inferred from the radial velocity 
data proved stable for inclinations down to i ~ 30°. At smaller inclinations, the masses of 
the super-earth planet s would become so large that the inner planet GJ 581e would have 
gotten ejected quickly. IBeust et al.l (120081 ) also computed the maximum Lyapunov exponents 
(MLE) and concluded that planet d is less chaotic than the inner planets. Within a wide 
range of initial parameters, the system is chaotic but long-term stable, which indicates that 



the regions of bounded chaos (jLaskarlll990l ) are extensive. 
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Fig. 2. — Orbital motion of the planet GJ 581d was integrated over 10 8 yr forward. This 
figure shows the first 10,000 years of evolution of this planet's eccentricity (red line) and 
semimajor axis (black line). 



Using the values listed in Table [T] as initial parameters, and assuming the orbits to be 
coplanar, we performed multiple simulations of the dynamical e yolution of the 5-body GJ 
581 system. The symplectic integrator HNBody, version 1.0.7, (IRauch fc Hamilton! 120021 ) 
was utilised with the symplectic option. We chose a time step of 5 • 10 -5 yr (0.018 days) 
and integrated the trajectories of all five bodies over 100 million years. As previous studies 
have found, the system appears to be quite stable. The eccentricities and semi-major axes 
show small periodic variations over the integration time. For GJ 581d, these variations over 
the first 10,000 years of integration are presented in Fig. |2j The eccentricity of this planet 
oscillates, seemingly forever, within a narrow range around 0.27. The semimajor axis, too, 
varies very little. 



In order to quantify the cha os in GJ 581, we employed the sibling simulation technique 
developed by iHayesI ( 120071 . 120081 ) to investigate the behaviour of the outer Solar system with 
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Fig. 3. — The distance between two sibling trajectories of the planet GJ 58 Id, with an 
initial perturbation factor in the semi-major axis of 10~ 14 , integrated over 10,000 yr. The 
separation grows exponentially just before 4000 yr, confirming the presence of chaos in the 
orbital motion of the system. 
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slightly different initial conditions. Two sibling trajectories are generated by perturbing the 
initial semi-major axis of the planet GJ 581b by a factor of 10 -14 . The distance between 
the unperturbed planets and their siblings is then computed as a function of time. An 
exponential divergence between the trajectories indicates that the system is chaotic. We 
find that chaos sets in within 5000 yr, with a characteristic Lyapunov time of only ~ 30 yr. 
The distance between siblings is shown for planet GJ 581d in Figure [3j In this case, the 
epoch of exponential divergence is close to 4000 yr, with a characteristic Lyapunov time of 
only 38 yr. The dominating planet b, on the other hand, displays two epochs of exponential 
divergence, with Lyapunov times of 18 and 277 yr. 

Long-term dynamical stability coupled with strong but bounded chaos is not uncommon 
among the currently known systems of multiple exoplanets. For example, a similar chaotic 
behavior was found for the dynamically robust system 55 Cnc, where a hidden commensu- 
rability may exist for the dominating 14-day planet. This makes the Solar system, whose 
chaos is much slower, to stand out. In the context of our study, the important conclusion is 
that the parameters of GJ 581d listed in Table [2] can be safely taken constant over extended 
intervals of time. 



4. Rotational evolution of GJ 581d. A combined effect of the triaxiality and 

tides 

-> (TBI) 

A planet is subject to the torque T due to the planet's permanent triaxiality, and 

-»(TIDE) 

to the torque T generated by tidal deformation. In our model, we shall assume that 
both these torques are exerted upon the planet by its host star only. We shall thus neglect 
the torques exerted upon the planet by its moons or by other planets. 



4.1. The equation of motion 

Consider a planet of the mean radius R and mass M p i anet and treat it as the pri- 
mary. The tide-raising perturber (the star of mass M sta r ) will be regarded as the secondary 
effectively orbiting the primary. 

The principal moments of inertia of the planet will be denoted as A, B, C , in assump- 
tion that A < B < C . The maximal moment of inertia will read as C = ^M p i anet R 2 , the 
numerical coefficient £ assuming the value of 2/5 in the homogeneous-sphere limit. 

Suppose that the planet is rotating about its major-inertia axis z , the one corresponding 
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to the maximal moment of inertia C . Let 9 be the sidereal angle of this rotation, reckoned 
from an arbitrary line fixed in inertial space (say, the line of apsides) to the axis x of the 
largest elongation, i.e., the axis corresponding to the minimal moment of inertia A of the 
planet. Rotation will then be described by the equation 

f-j- (TRI) (TIDE) (TRI) —(TIDE) 
Q _ J_z ' ' z _ l_z ~r I z 

a i let 

with the subscript z serving to denote the polar components of the torques. High-accuracy 
integration of equation (J3J) requires a step much smaller than the orbital period. 



4.2. The triaxiality-caused torque 

We approximate the polar component of the torque with its quadrupole part given by @ 



= 3 -(B-A)^sin2^ (5a) 



« - - (B- A) n 2 ^ sin2(0-/) . (5b) 

Here r denotes the distance between the centres of mass of the two bodies, while M star 
stands for the mass of the star (which, in our setting, is playing the role of perturbing 
secondary effectively orbiting the tidally deformed primary). The mean motion is given by 
n = a/ (M p i anet + M star ) I a 3 w \/M star /a 3 , in understanding that the star is much more 
massive than the planet. 

Our notations for angles are given on Figure HI The angle ift renders separation be- 
tween the planetocentric direction towards the star and the principal axis x of the maximal 
elongation (the minimal-inertia axis). It is equal to the difference between the angles / and 
9 made by these two directions on an arbitrary line fixed in inertial space. If this fiducial 
line is chosen to be parallel to the line of apsides connecting the star with the empty focus, 
then / will be the planet's true anomaly as seen from the star, while 9 will be the sidereal 
angle (which we agreed to reckon from the line of apsides.). 



Employing f!5bj) in practical computations, it is advantageous to approximate the func- 
tions (a/r) 3 sin 2 f and (a/r) 3 cos 2f with truncated series of the Hansen's coefficients (e.g., 



Branhamlll990l ). which in their turn can be approximated with series of the Bessel functions 



of the first kind. 



See, e.g.. lDanbvl (|1962() 
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Fig. 4. — The horizontal line is that of apsides, so / is the true anomaly, 9 is the sidereal 
angle of the planet, while if) = f — 9 . The principal axes x and y of the planet correspond 
to the minimal and middle moments of inertia, appropriately. 
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4.3. The tidal torque 

Due to a several-orders-of-magnitude difference in the intensities of internal friction in 
a terrestrial planet and a star, we shall take into account only the tides exerted by the star 
on the planet, and shall neglect the tides on the star. 

The technique of calculation of tidal torques had until recently remained in a somewhat 
embryonic state requiring much of correction. On the one hand, it had long been habitual and 
common to combine the expression for the tidal torque with unphysical rheological models. 
As a result, much of the hitherto obtained results on the spin history of terrestrial planets 



[e.g. 



Celletti et al.ll2007l ; iHeller et al.ll201ll ) have to be reexamined, because they were based 
on very ad hoc rheologies incompatible with the behaviour of realistic solids. On the other 
hand, in many publications the description of the tidal torque was marred by a mathematical 
mistake which proliferated through many papers and ended up in textbooks. Expl a natio n 



and correction of that long-standing oversight is presented in lEfroimsky fc Makarovl (120121 ). 



We shall limit our treatment to the simpler case where the planet is not too close to the 
star (R/a 1), the obliquity of the planet is small (i ~ 0), and its eccentricity e is not very 
large. As demonstrated in the Appendix, under these assumptions the polar component of 
the secular part of the tidal torque can be approximated by 

(7 - (TIDE)) = 

X ' 1=2 

a 7 

- G M s ] ar R 5 a" 6 G 2 (e) sin | e,^) | Sgn ( u 2Wq ) + 0(e 8 e) + 0(z 2 e) , (6) 

where R denotes the radius of the planet, 9 stands for its rotation rate, while a , e , and 
% are the semimajor axis, the eccentricity, and the obliquity of the perturber (the star) as 
seen from the planet. The notation G\ vq [e) stands for the so-called eccentricity polynomials 

(-1-1), (l-2p) 

which coincide with the Hansen polynomials +g \ (e) . 



For the first time, this expression appeared, without proof, in the work by lGoldreich fc Peale 



( 119661 ) who summed over all integer values of q . For the reasons explained in Appendix A 
below, we limit the summation to the span from q = —1 through q = 7. This truncation 
effectively limits the Taylor expansions of eccentricity-dependent functions to terms of the 
orders up to e 7 , inclusive. The resulting relative precision of our calculations is approxi- 
mately 0.1%, because the largest and the smallest terms of the so-truncated sum contain 
G 200 (0.27) = 0.82202 and G 20 7(0.27) = 0.10206 , respectively; while all the eccentricity poly- 
nomials G2o q (e) outside the range of summation (i.e., for q outside of the range — 1, . . . , 7) 
assume much smaller values. 
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In equation ([6]), both the dynamical Love number k 2 and the phase lag e 2 are functions 
of the tidal Fourier mode which in this, simplified case reads as 

^220? = (2 + q)n - 26 , (7) 

a more general expression for ui mpq being given in the Appendix. As ever, n denotes the 
mean motion. 

The dynamical Love numbers k 2 (u 220 q) are positive definite, while the sign of each lag 
e ^( u 22oq) coincides with that of the tidal mode oo 220 q , as explained in Appendix A. Therefore 
each factor k 2 (u 220 g) sin e 2 (cJ 220 g) can be written as k 2 (u 220 q) sin | e 2 (w 220 g) I Sgn ( UJ 22oq) ■ 

Through (J7J), each such factor becomes a function of the planetary spin rate 9 : 
H^wq) sin I e 2 (u 220q ) | Sgn ( u 220q ) 

= k 2 (2(n - 9) + qn) sin | e 2 ( 2(n - 6) + qn ) | Sgn ( 2(n - 9) + qn ) . (8) 

Accordingly, the entire sum ([6]) can be regarded as a function of 9 . The mean motion and 
eccentricity act as parameters, their evolution being much slower than that of 9 . 

In each term of series flSD, the factor k 2 sine 2 , expressed as a function of 9, has the 
shape of a kink. In Figure the dotted line □ illustrates the behaviour of the factor 

^2202) sine2 K 2 o 2 ) = fc 2(4n - 20) sine 2 (4r2 - 2 9) 

= k 2 (An-29) sin|e 2 (4n - 29) \ Sgn(4n - 29) . (9) 

The origin of the kink shape is explained in Appendix B, with references provided. 
From the physical viewpoint, the emergence of this shape is natural, in that each term 
should transcend zero and change its sign smoothly when the spin rate goes through the 
appropriate spin-orbit resonance. 

However, in ([6]) the kink-shaped factors are accompanied with different multipliers 
G 2 (e) . So the sum as a function of 9 , will be a superposition of many kinks having 

different magnitudes and centered at different resonances (eight kinks, if the sum over q 



5 In a hypothetical case of a planet despinning at a constant rate through a resonance, the appropriate 
tidal mode becomes linear in time. Then the tidal torque assumes a similar kink shape, as a function of 
time. This situation is considered in Ferraz-Mello (2012, Figure 7b). Any physically reasonable rheology 
must lead to this or similar type of tidal torque behaviour in the vicinity of a resonance. 
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goes from -1 through 7). The resulting curve will cross the horizontal axis in points close to 
the resonances 9 = n (1+ q/2) , but not exactly in these resonances — see the solid line in 
Figure EJ 



Goldreich &; Peald (119681 ) were first to point out that in the vicinity of a particular 
resonance q' , i.e., for 9 being close to n (1 + q'/2) , the right-hand side of ([6]) can be 
conveniently decomposed into two parts. One part is the q = q' term, a kink-looking odd 
function vanishing at 9 = n (1 + q'/2) . Another part, called bias, is comprised by the rest 
of the sum. So the bias is the contribution of all the q ^ q' modes into the values assumed 
by the torque in the vicinity of the q = q' resonance. For not too large eccentricities, the 
bias is usually negative in value. The bias is a very slowly changing function, which can, to 
a good approximation, be treated as constant. 

The q = q' term by itself is an odd function, and it goes through nil at exactly 
9 = n (1 + q'/2) . However, the bias displaces the location of zeroes. In Figure |5j the 
torque (depicted with a solid line) is defined mostly by the term with q = 2 (rendered by 
the dotted line). However, the curve is shifted down due to the bias which is defined mainly 
by the right slope of the q = 1 kink located to the left. As the right slope of the q = 1 kink is 
negative, the q = 2 kink in Figure |5] is shifted slightly down, and the zero is located slightly 
to the left of the point 9 = 2n . The crossing point in this case is at 9/n = 1.999976481 . 
Such a minuscule shift of the equilibrium away from the resonance frequency does not have 
any practical consequences, because the net nonzero tid al torque is compensated by a tiny 



secular triaxial torque, as discussed in more detail in ( IMakarov fc Efroimskyl 120121 ) . The 
shifts of the tidal-equilibrium frequencies at resonances are larger for bodies with a lower 
Maxwell time, e.g., for bodies whose mantles contain a large fraction of partial melt. 



5. Probabilities of capture into resonances. General facts 



When the planet's spin rate approaches a resonance 29 = {2 + q)n , with an integer q, the 
planet can be captured, or can traverse the resonance, depending on the specific trajectory 
in the phase space. A method for estimation of capture probabilities was developed by 
Goldreich fc Peald ( 119661 ). for two simpli fied models of the tidal tor que. A clearer explanation 
of this theory was presented in the later I Goldreich fc Peald (119681 ) to which we shall refer. 



In one model considered in these papers, the torque was assumed to be linear in the 
tidal frequency. To be more precise, in the Fourier expansion of the torque over tidal modes, 
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Fig. 5. — Angular acceleration of the planet GJ 581d caused by the secular tidal torque (EJ), 
in the vicinity of the 2:1 spin-orbit resonance. The dotted kink-shaped curve depicts the 
term q = q' , which is an odd function when centered at 6/n = 1 + q'/2 . In the paper, 
we denote this term with W{^) . The solid line furnishes the overall torque calculated as a 
sum of the q = q' kink and the bias. The bias, denoted in the text with V , is comprised 
by the terms of <Q, with q ^ q' . 
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each torque co mponent was set to be p roportional to the appropriate tidal mode § 
formula (19) in lGoldreich fc Pealei (I19681 ). 



- see 



Another model addressed in the said two works was the constant-torque one. Specifi- 
cally, in the Fourier expansion of the torque over tidal modes, each torque component was 
set to be a constan t mult iplied by the sign of the corresponding mode - see equation (29) in 
Goldreich fc Peale 



In their treatment of both models, iGoldreich fc Peald (119681 ) took into account only 
the secular, orbit-averaged, component of the tidal torque, and ignored the existence of an 
oscillating component. Below we shall test the validity of this approximation. 

The pivotal ideas and formulae of the capture theory are explained in short in Appendix 
C below. Here we shall employ some o f those formulae, th o ugh w ith an important difference. 
Instead of the toy models employed in lGoldreich fc Pealei (119681 ). we shall rely on a realistic 
rheology of solids. As we shall see, capture probabilities are sensitive to changes in (at least 
some of the) rheological parameters. For example, the probabilities turn out to be more 
sensitive to the mantle's Maxwell time than to the planet's triaxiality. 



The principal result of the analysis carried out by IGoldreich fc Peald (119681 ) is their 
estimate for the probability of capture into an arbitrary resonance q = q' , i.e., into a 
steady rotation at the rate of 9 = n (1 + q ' /2) . The probability is given by 

2 



capt 



l + 27rV/f: n W(<y)d 1 
where the new, auxiliary variable 7 is defined through 



7 



2 9 - (2 + q')M 



its time derivative thus being the negative double of the tidal mode <^ 2 2oq' '■ 



7 



- (2 + q ')n + 29 



b~>220q> 



(10) 



(11) 



(12) 



This mode vanishes in the q = q' resonance, so we may say that this resonance corresponds 
to 7 = . The odd function W{^) entering expression (fTOl) is called kink and is given by 



W{fi) 



K G 22oq> k ^ sin l e 2(7)l Sgn(7) 



(13) 



6 The term tidal mode is more appropriate than frequency, because a Fourier mode can assume either 
sign, while the physical frequencies are the modes' absolute values and thus are positive definite. 



7 As demonstrated in (|Efroimskv fc Makarovll2012l) . this model contains an inherent mathematical con- 
tradiction . 
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K bein g a positive cons tant. Be mindful that our defini tion of 7 is in agre e ment with that 
used by iMakarovl (120121 ) and is twice that introduced in lGoldreich fc Peald (119661 Il968l ). 

Entering formula (1101) is also a function 



V 



J 220q) S111 I e 2l W 220gJ 



Sgn (cu, 



K ^ G 22oq hdq - q')n) sin|e 2 ((g - q')n) | Sgn (q - q') 



(14) 



called bias. The physical meaning of V and ^(7) is explained in Appendix C. The 
Appendix also explains the way how in the above integral the dependence of 7 upon 7 
should be set, so that the integral could be evaluated. 

Important to us is the fact that the above expressions for V and ^(7) contain the 
mode-dependent factor /c 2 (w 2 2og) sm e2 ( a; 22og) • The functional form of its mode- dependence 
is defined by the self-gravitation and rheology of the planet (jEfroimskyl l2012al ). This ob- 
servati on helps us to mark the point at which our analysis will diverge from that carried 
out by iGoldreich fc Peald ( 119681 ) : our choice of a realistic rheology will yield for /c 2 sin e 2 a 
mode-dependence different from both cases addressed in Ibid. 
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Fig. 6. — Probabilities of capture of GJ581d at 3:2, 2:1, 5:2, and 3:1 spin-orbit resonances 
as functions of eccentricity, for (B — A)/C = 5 • 10~ 5 . Left: a warmer planet, with a lower 
viscosity and tm = 50 yr. Right: a colder planet, with a higher viscosity and tm = 500 yr. 



8 Writing the right-hand side of (jl4|) , we used the fact that we are in the vicinity of the q = q' 
see formula (|39|) in Appendix C. 



resonance 
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6. Probabilities of capture into resonances. The case of a realistic rheology 



Recently, iMakarovl (120121 ) generalised the treatment by iGoldreich fc Peald (Il968[ ) to the 
case of a more physical tidal torque, i.e., of a torque whose terms bear a mode- dependence 
stemming from the properties of realistic solids. The rheology of silicates and ices is well 
described by the Andrade model, insofar as the tidal frequency is sufficiently high and the 
inner friction is dominated by anelasticity. At lower frequencies, the friction is dominated 
by viscosity, and t he behaviour of the material becomes closer to that of the Maxwell body 
( lEfroimskyi l2012aU bl) . The resulting mode-dependencies of the terms of the torque are pre- 
sented below in Appendix B. These mode- dependencie s are more complicated than in the 
two toy models considered by IGoldreich fc Peald (119681 ) . 



In our computations, we ignore the oscillatory components of the tidal torque, and 
consider only the secular part which can be obtained by averaging over one orbital pe- 
riod — a convenie n t sim plification is to be justified below. Thus we begin in the spirit of 
Goldreich fc Peald (Il968[ ). i.e., use formula (fTUj) . the f unctions W an d V being furnished 
by formulae ( TT3|) and ( !T4|) . correspondingly. Following IMakarovl (120121 ). we then insert into 
those formulae the realistic mode- dependence of &2 sine2 written down and explained in 
Appendix B. Finally, we perform a brute force numerical check of whether the neglect of the 
oscillating part is legitimate. 

The resulting capture probabilities as functions of eccentricity are presented in Fig. [6] 
for two values of the Maxwell time: r M = 50 yr (the left plot) and r u = 500 yr (the right 
plot). While 500 yr is the current Maxwell time of the Earth's mantle, the value of 50 yr is 
deemed to represent a slightly % warmer planet, one with a lower viscosity. 

From the plots, we see that probabilities of capture are sensitive to the Maxwell time. 
Warmer planets with lower tm , have more chances of being captured at spin-orbit reso- 



9 We say slightly , because the viscosity rj depends upon the temperature T through the Arrhcnius formula 
r] oc cxp(A* / RT) , while the rigidity [i depends upon T slower, unless we are close to the melting point 
(the latter caveat being hardly relevant, as the lithostatic pressure prevents the mantle from inciting even 
though some partial melt may be present). So the Maxwell time t m = rj/fi depends on the temperature 
about exponentially. A relatively small variation of the temperature renders the following relative changes 

At At? AT 1 4* 

of the viscosity and the Maxwell time: — _ M ss ^ w . Here the activation energy may be 

M '/ 1 It 

estimated, for olivines and silicate perovskites, as A* w 6 x 10 5 J/mol. Consider an Earth-like planet with 
a silicate mantle of a mean temperature T = 2300 K. A decrease of the viscosity and the Maxwell time by 
9/10 would correspond to an increase of the temperature by AT = 66 K. For the Earth, it would imply a 
less than 1 Gyrs step back to the Neoproteozoic Era, one marked with appearance of the first multicelled 
organisms. Therefore both values, 500 yr and 50 yr, serve as legitimate estimates for the Maxwell time of 
an Earth-like, potentially habitable planet. 
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nances, as they loose their spin angular momentum. In particular, for tm = 50 yr and 
(B — A)/C = 5 ■ 10 -5 , the capture probabilities are: 1 in resonance 3:2, 0.897 in 2:1, 0.383 
in 5:2, 0.133 in 3:1, and 0.040 in 7:2. Therefore, if GJ 581d has had enough time to de-spin 
into a state of spin-orbit equilibrium at the current value of eccentricity, then the probabil- 
ities of the planet being now in resonances are the following: 0.05 in resonance 3:2, 0.46 in 
2:1, 0.32 in 5:2, 0.13 in 3:1, and 0.04 in 7:2. 

For tm = 500 yr and (B — A)/C = 5- 10~ 5 , the capture probabilities are: 1 at resonance 
3:2, 0.568 in 2:1, 0.188 in 5:2, 0.054 in 3:1, and 0.014 in 7:2. Then, with a probability of 0.35, 
the planet is currently entrapped in resonance 3:2, 0.43 in 2:1, 0.18 in 5:2, 0.05 in 3:1, and 
0.01 in 7:2. The 2:1 resonance is the most likely state of the planet in both cases. However, 
for tm = 50 yr the second likeliest state is 5:2, while for tm = 500 yr it is 3:2. 

The probabilities of capture depend on the degree of triaxiality through the parameter 
(B— A) /Cm the equation of separatrix (1371) . This equation tells us that larger values of (B — 
A)/C entail greater libration amplitudes of 7 . The kink W^) , which is the antisymmetric 
part of the secular torque, becomes smaller with 7 growing outside the resonance, see Figure 
[5J Therefore, the integral f_ ^(7) d'y in equation [42] is expected to become smaller for 
larger (B — A) jC . Accordingly, the capture probability is to become smaller for larger 
(B — A)jC . This is confirmed by computations of capture probabilities for tm = 50 yr, 
(B — A)/C = 2 • 10~ 4 and the other parameters as given in Table [2j The capture probabilities 
are: 1 in the resonance 3:2, 0.776 in 2:1, 0.301 in 5:2, 0.097 in 3:1, and 0.028 in 7:2. Each of 
these probabilities is slightly smaller than its counterpart taken for the smaller triaxiality 
(B — A)/C = 5 • 10~ 5 . At the same time, the relatively small differences between the capture 
probabilities into the same resonance, for different values of (B — A) jC , imply that the 
spin-orbit resonances are more sensitive to tm than to [B — A)/C . 

As a way of spot-check verification of these theoretical results, we conducted brute-force 
simulations of the GJ 58 Id system. The equation of motion, incorporating both the tidal 
and triaxial torques acting on the planet, was numerically integrated 40 times for tm = 50 
yr, (B — A)jC = 5 ■ 10 -5 , e = 0.27, the initial spin rate #(0) = 2.51 n, the initial mean 
anomaly .M(O) = 0, and the initial sidereal angle 0(0) = ni/40, with i = 0,1,..., 39. 
This method of estimation tacitly assumes that the initial sidereal angle at a fixed rate of 
rotation can take any values with equal probability (which is less than obvious). With this 
assumption accepted, simulations spanning 7000 yr, with a step of 1.5 • 10~ 3 yr, resulted in 



10 For terrestrial planets, as well as for large solid moons, our choice of the values 5-10 5 and 2 • 10 4 
for (B - A) jC is likely to be realistic. Recall that (B - A) jC is equal to 2.2 ■ 10 -5 for the Earth (Wen-Bin 
Shcn et al. 2011, Liu & Chao 1991), to 6.9 • 10~ 4 for Mars (Edvardsson et al. 2002), and to 2.3 • 10" 4 for 
the Moon (Williams et al. 1996). 
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14 captures into the 5:2 resonance and 26 passages. The estimated probability of capture 
is thus 14/40 = 0.35 , which is surprisingly close to the theoretical estimate 0.383 . To 
make the numerical and theoretical estimations consistent, the former included only the 
secular part of the torque (P). This points at one of the weaknesses of the semi-analytical 
derivation of probabilities. The oscillatory terms of the tidal torque are ignored altogether. 
Furthermore, even though the defor mity of the planet ( B — A )/C enters the computation 



of capture probabilities, performed by lGoldreich fc Peald (119681 ) . the cyclic variations of the 



spin rate caused by the triaxial torque are not involved in any way. The smooth sinusoidal 
separatrix trajectories defined by equation (l3Tj) are in fact superimposed with a jitter caused 
by the harmonics of time-dependent terms of the torque. 

Figure [7J illustrates the role of the oscillatory torques which are ignored in the theoretical 
calculation of capture probability. The graph renders the behaviour of the quantity 7/2, 
which is proportional to the energy of rotation, as a function of 7 for two full libration 
cycles, one going forward toward the point of resonance ( 7 > ) and the other going back 
toward 7 = beyond the point of resonance (7 < 0). Due to the dissipation of energy by 
the tides, the latter curve is systematically lower than the former one, but the difference is 
so small that it cannot be seen on the graph. Besides, the oscillatory terms of the net torque 
make the curves jittery. Although the amplitude of the jitter diminishes in the vicinity of 
the resonance, it may app ear to be significant for the outcome of this process. According to 



Goldreich fc Peald ( 119681 ). capture occurs when the first post-resonance minimum of 7 2 /2 
reaches 0. If the jitter is superimposed with the smooth separatrix trajectory, the chances 
to bump into seem greater than without jitter. The inset in Fig. [7J shows in much greater 
detail this important segment of the post-resonance libration, in the axes lg(7/2) versus 7 . 
The rotational energy comes very close to at the lower extent of oscillations, never quite 
reaching it. And indeed, in this simulation, the planet traversed the 2:1 resonance. 

To get a better understanding of the influence of oscillating terms of tidal torque on the 
chances of the planet to be captured, we repeated the 40 high-accuracy integrations described 
in the previous paragraph, at the same initial parameters, though this time including the 
entire set of periodic terms. We found that the results for individual simulations often 
changed, i.e., what resulted in capture for a given set of parameters became a passage, 
and vice versa. Surprisingly, we recorded 14 captures out of 40 trials, yielding the same 
probability of 0.35. From this small-scale experiment, we see that, while the outcome of 



11 Recall that the common expression ([6]) for the polar component of the tidal torque renders only its 
secular part. It is for this reason, that in this formula we use the angular brackets denoting an average over 
an orbital period. A full expression for the torque includes also an oscillating pa rt which averages out over 
an orbital period, but which may nonetherless play a role in the capture process (|Efroimskvll2012bf ). 
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Fig. 7. — Two librations immediately before and after the point of 2:1 resonance. The 
resonance is traversed in this case, because the descending branch of the second libration 
never reaches 7 = 0. The inset shows in more detail a fragment of the crucially important 
minimum of the post-resonance libration in logarithmic scale. 
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a particular integration may change because of the jitter of the tidal torque, the overall 
probability of capture is unlikely to depend on it. 



7. Conclusions 

The growing number of detected systems of multiple exoplanets and the impressive 
quality of observational data collected for them allows astronomers to perform analysis of 
probable dynamical states and evolution of these remote worlds at a level of detail unthink- 
able just several years ago. Still, this analysis is riddled with difficulties and uncertainties. 
The planets of G J 581 present a challenge for both observational practice and interpretational 
theory. The story of two fictitious planets tells a lesson about the hazards of combining, with- 
out proper caution, the data from two instruments with their own sets of systematic errors. 
It also calls for a certain unification of the planet detection techniques or, at least, for an 
easily accessible and well-tested standard detection algorithm available as a web application. 
It is fine to apply a variety of different detection methods to the same data, but the standard 
algorithm should always be checked, and discrepancies, if any, should be investigated and 
reported. iV-body integration of detected systems should also be a norm reducing the prob- 
ability of blunder. The system of G J 581 proves to be remarkably stable, with the four bona 
fide planets remaining on their orbits despite the strong evidence of chaos. The characteristic 
Lyapunov times are very short compared to the dynamical lifetime of the system. 

We have explored the rotation history of the planet GJ 581d assumed to have compo- 
sition alike to that of the terrestrial planets of the Solar system. 

Contrary to the previous publications on the subject, which were based on ad hoc, 
simplistic tidal theories, we find that this planet cannot be captured into a synchronous or 
pseudo-synchronous rotation, if it began its evolution from faster, prograde initial spin rates. 
Instead, for a plausible range of parameters, the most likely state of the planet is the 2:1 
spin-orbit resonance. In this state, the day on GJ 581d should be 67 Earth's days long, 
which bides for an inhospitable environment, though the potential habitability of this planet 
cannot be ruled out just on climatic considerations. The case of 2:1 spin-o rbit resonance was 



consid ered in the simulations of a hypothetical atmosphere on GJ 581d by I Wordsworth et al. 



(120111 ) whose modeling confirmed the possibility of liquid water being present on the surface, 



under some favorite conditions. 



The next likeliest equilibrium states are the 3:2 or 5:2 resonances, depending on the 
temperature and viscosity of the mantle (much less on the planet's triaxiality) . 



At the same time, in the event that the initial rotation of the planet was retrograde, 
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the most probable final state is synchronous rotation. 
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Appendix A. The tidal torque 

As well known, the tidally generated amendment U to the potential of the perturbed 
planet can be presented in the form of a Fourier series over the tidal modes 

Lo lmpq = (l - 2p) Co + (I - 2p + q) M + m(tt - 9) « (I -2p + q)n - m0 , (15) 

where 6 denotes the rotation rate of the tidally perturbed primary (the planet), while u , 
Q , n , and M. are the periapse, the node, the mean motion, and the mean anomaly of the 
perturbing secondary (the star) as seen from the primary. While the tidal modes oo lmpq can 
be of either sign, the actual physical forcing frequencies Xi m pq a ^ which the strain and stress 
oscillate in the perturbed body are positive-definite: 

Xlmpq = | UJlmpg | « | (I - 2p + q) 71 - TU 6 | . (16) 



The full ser i es for the tidally generated amendment to the planet's potential was written 
down by iKaulal (|1964[ ). though its partial sum was known yet to Sir Charles Darwin (1879). 
For this reason, we often term this Fourier series as the Darwin-Kaula expansion. We apply 
this term also to the ensuing series for the tidal torque acting on the perturbed planet. 

A detai led derivation of th e Fourier expansion of the polar component of the torque can 
be found in lEfroimskyl (I2012bf) . It turns out that the torque contains both a secular and a 
rapidly oscillating part] 12 ! the secular part being given by 



R 



21 + 1 



(I 



m) 



21 + 2 XI ( 

1=2 CI m=0 V >' 



I 

p=0 



(^Impq) Sin €i {UJlmpg) , (17) 



12 The oscillating part averages to nil and is not expected to reshuffle much the probabilities of capture 
(though it can effect the fate of each particular trajectory). It can however influence the process of damping 
of free librations. 
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where G is Newton's gravity constant, a, z, e denote the semimajor axis, inclination, and ec- 
centricity, while the angular brackets (...) signify orbital averaging. The standard notations 
Fimp(i) and Gi pq (e) are used for the inclination functions and the eccentricity polynomials. 
The Love numbers k [ and the phase lags e l are functions of the Fourier tidal modes u) lmpq 
given by (fIBl) . 



As explained, e.g., by lEfroimsky &; Makarovl (120121 ). in the Darwin-Kaula theory of tides 



the phase lags emerge as the products of the modes ooi mpq by the appropriate time lags: 

e l{Ulmpq) = Uimpg Ati((J[ mpq ) , (18) 

where the lags are written down not as £i mpq and Ati mpq but as ei(ui mpq ) and Ati(ui mpq ) . 
The same concerns the notation for the Love numbers. This is done in order to emphasise that 
for a homogeneous near-spherical body the functional forms of the lags and Love numbers (as 
functions of the Fourier mode) are defined solely by the integer number I (called degree) . Ef| 

For causality reasons, the time lags At z (o; /mp(? ) are positive-definite. Therefore, f[TB"j) 
may be rewritten as 

tli^lmpq) = Xlmpq Atl(oJl mpq ) Sgn ( Ul mpq 

) , (19) 

Ximpq being the physical forcing frequency (TT6]) . As a result of this, the entire expression for 
the polar component of the torque can be written down as 

(7 -(TI DE ) ) = 

_^2l + l i i 

2 GM s 2 tar ~2T+2~ (I T ra,, 771 F l 2 mpi % ) Y G U e ) k ^ U lmpq) I ^Impq) \ Sgn ( U} lmpq 

1=2 & m=0 ^ *' p=0 q=-oo 



13 



In his cornerstone work, iKaulal (|1964 ) employe d the somewhat incons i stent notations ki and e 



' ' ' ■ - I 1 I 1 Impq 

which were later borrowed by other authors, e.g., by Efroimskv fc Williams! ( 20091 ) who also used a similar 



notation Ati mpq for the time lag. In our later works ( Efroimsky||2012al rJ). we chose to switch to e; and Ati 



since the forms of the functional dependencies of the lags on the modes are defined by the degree 

I only. The lags' (and Love numbers') dependence on the other three integers, m, p, q , originates only 
through the dependence of the argument w impg upon these integers. This is why the notations fc;(w impg ) , 
e i(ui mpq )> ^li^impq) are more adequate than fc ( , e lmpq and At lmpq . 

All said applies only to near-spherical homogeneous celestial bodies. For nonspherical bodies, the 
situation becomes more involved, as the functional form of the lags and Love numbers acquires dependence 
on all the four integers. In this case, we should write k lmpq , e ; , and At lmpq . Fortunately, for a slightly 
non-spherical body, the Love numbers and lags differ from the Love numbers and lags of the spherical 
reference body by terms of the order of the flattening, so a small non-sphericity can be neglected. 
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When the planet is not too close to the star (R/ a <C 1), it is possible to neglect the terms 
with / > 2. For small obliquities (i ~ 0), it is also possible to leave only i— independent 
terms (the next-order terms being quadratic in i). Finally, for not too large eccentricities 
( e 1 ), it is reasonable to take into account only the terms up to e 7 , inclusive. This would 
imply summation over the values of q running from — 7 through 7 . In reality, however, the 
term with q = — 2 turns out to vanish identically, while the terms with q = — 7 through 
q = — 3 are accompanied with extremely small numerical factors and can thus be dropped. 
This way, the polar component of the tidal torque gets approximated with 

(T 2 (TIDE) > = 

x 1 1=2 

- G M s 2 tar R 5 a-* G 2 {e) Huj^) sin\ e,^) | Sgn ( u 2Wq ) + 0(e 8 e) + 0(z 2 e) , (21) 



Histor ically, this expression first appeared, without proof, in the paper by lGoldreich fc Peak 



(119661) who summe d over all integer values of q . A schematic proof was later offered by 



Dobrovolskisl (120071 ). 



The shape of the factors ki sine; as functions of the mode 0J lmpq is defined by the 
size and mass of the body and by its rheology. A rheological model is a constitutive equa- 
tion interconnecting the strain and stress. Within linear rheologies, such equations can be 
rewritten in the frequency domain where each harmonic of the strain gets expressed alge- 
braically t h rough the appropriate harmonic of the stress. Using the techniques explained in 
Efroimsky (l2012al lbl). those algebraic relations can be employed to derive the shape of the 



functions ki(u lm ) smei(u lmpq ) standing in the terms of the Darwin-Kaula expansion of 
the tidal torque. 

We rely on the Andrade rhe ological model which is known to work for the Earth's 



mantle (lEfroimsky fc Laineyl 120071 ) and which may therefore be applicable to the mantles of 
other terrestrial planets. The universality of the Andrade model lies also in the fact that it 
can be rewritten i n a manner permitting a switch to the Maxwell model at low frequencies 



jEfroimskv! [2012al Jbll The necessity for this switch is dictated by the fact that different 



physical mechanisms dominate friction over different frequency bands. As demonstrated in 
Ibid., employment of this combined model (Andrade at higher frequencies, and Maxwell at 
lower frequencies) renders for ki sin e; a kink-shaped dependence upon the Fourier mode - 
as the dotted line on Figure [5j 

In practical calculations, it is convenient to insert ffT5]) into ki(u; l ) sin e\ (ou lmpq ) 
and thus to obtain the dependencies of k\ sin q upon the spin rate 9 , with n treated as 
a constant or slow- varying parameter. The dependence of fc; sine; upon 9 will have the 
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shape of a kink too. It will be a function varying slowly everywhere except in the vicinity 
of the spin-orbit resonances 9 = - — ^ - n . The dotted curve in Figure \5\ depicts the 
dependence of the factor 

H^ 2202 ) sine 2(w 2202 ) = k 2 {An - 2 9) sine 2 (4n - 20) 

= k 2 (An - 29) sin|e 2 (4n - 29) | Sgn(4n-20) (22) 

in the vicinity of the 2:1 spin-orbit resonance. This factor shows up in the 2202 term 
of the torque. In the 1:2 spin-orbital resonance, 9 transcends the value of 2n , so the 
tidal mode oj 22 $ 2 = 2{2n — 9) goes through nil and changes its sign. In Figure the 
factor k 2 (u 2202 ) sin e 2 (w 2202 ) does the same: as the tidal mode u 2202 approaches zero (or, 
equivalently, as 9 goes through 2n), the said factor smoothly goes through zero and changes 
its sign. This makes the considered term of the tidal torque change its sign smoothly when 
the synchronous orbit gets crossed. Note that the rapid but smooth changes of tidal torque 
take place within a very narrow interval of 9 . The widely used assumption that the torque 
component is linear in 9 (and therefore linear in the tidal mode) is not justified by this 
model, except in an extremely small range of spin rates around the point of resonance. 

Similarly, for any set of integers Impq , the factor 

H^lmpq) Sine *Kn W ) = H^lmpq) Sil1 I ^ ("impq ) I S § n ( "impq ) = 

ki( (/ — 2p + q) n — m9) sin | ei ( (I — 2p + q) n — m 9 ) \ Sgn ( (/ — 2p + q) n — m 6 1 ) , (23) 

depicted as function of 9 will demonstrate behaviour similar to that of (1221) in Figure [5j it 
will smoothly go through nil and will change its sign as the Impq commensurability gets 
transcended, i.e., as the tidal mode u) lrnpg = {I — 2p + q)n — m9 goes through nil. 

As ensues from (120]) and (121!) , the Impq term of the torque is decelerating for u lmpq > 
and is accelerating for u lmpq < . This motivates us to use the sign convention as in Figure 
the Impq term of the torque is agreed to be negative on the right of the Impq resonance 
and positive on its left. 

The overall tidal torque ff20|) . or its approximation (|2~T1) . taken as a function of the the 
spin rate 9 , will look as a superposition of kinks. In other words, if we sum up all the terms 
in (|2"U1) or (|2"TT) . and depict the sum against 9 , we shall get an overall curve containing a kink 

near each Impq resonance, i.e., near the points 9 = ^ — - n . We write near, because 



14 



This sign convention is opposite to the one we used in lEfroimskvl ([2012 
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these kinks will not go through zero at the said points exactly, but will be slightly displaced. 
This will happen because the higher-resonance kinks will be residing on the slopes of lower- 
resonance kinks. An example of this is furnished by Figure where the kink originating from 
the Impq = 2202 term is depicted with the dotted line. The total torque, i.e., the sum of 
the kink with the bias, is given by the solid line. We see that the presence of the bias shifts 
the kink slightly downward. So the total torque (the solid line) crosses the 9/n axis in a 
point located slightly to the left of the resonance 9/n = 2 . 



Appendix B. How rheology enters the play 



A laborious calculation (presented in lEfroimskyll2012al Jbl) demonstrates that the factors 
ki(u lmpq ) smei(u lmpq ) can be expressed via the real and imaginary parts of the complex 
compliance of the mantle material, and the mass and radius of the planet. This way, the 
mode-dependence of these factors is defined by both the rheology and self-gravitation of the 
planet. The factors come out to be odd functions, which is very natural, since each such 
factor (and the term of the torque, which contains this factor) must change its sign when 
the appropriate commensurability is transcended. 

Being odd, the factors can then be written down as k(u lmpq ) sin | ei(u lmpq ) | Sgn ( u lrnpq ) , 
the product ki(u lmpq ) sin | ci{uj 1 q ) | being an even function of the tidal mode. In other 
words, this product may be regarded as a function not of the tidal mode 0J lmpq but of its 
absolute value Xi mp q = I u im Pq I > wmcn is the physical forcing frequency of tidal oscillations 
in the mantle. All in all, we have: 



k M mpq ) shie,(u; ZmOT ) = h{u lmpq ) sin | e t (u lmpq ) | Sgn (u 



Impq 



= HXimpq) sin|e/(x impg ) I Sgn(w, mM ) . (24) 
The calculation in Ibid, furnishes the following frequency dependence for these factors: 



,, , . , , 3 -A, Jim [J(y) J 

M<0 ™«(<0 = W ^ Y) {Ke[m +Ajf - (Im[mr s ^«> • < 25 > 

where x is a shortened notation for the frequency X [mpq , while coefficients Ai are given by 



(2 I 2 + Al + 3)/i 3 (2 I 2 + 4/ + 3)/i 
IgpR : : 4 IwGp 2 R 2 



A i = rm? = — a i z n ,2 p2 — . ( 26 ) 



with R , p , a , and g being the radius, mean density, unrelaxed rigidity, and surface gravity 
of the planet, while G being the Newton gravitational constant. 
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The functions TZe[J(x)} and Im[J(x)] are the real and imaginary parts of the complex 
compliance J(x) of the mantle. These are rendered by the formulae 



Ke[J( X )] = J + J(xr A y a cos ( ^ ) r(a + 1) (27) 



2 



and 



lm[J( X )} = - J (xrj- 1 - J( X r A r a sin(y) I> + 1) , (28) 

J = 1/fi being the unrelaxed compliance of the mantle, and a being a numerical parameter 
assuming values of about 0.3 for solid silicates and about 0.14 — 0.2 for partial melts. In 
our computations, we used a = 0.2 . 

Among the rheological parameters entering ( 1271 - 128]) is the viscoelastic (Maxwell) time 
r M , which is the ratio of the mantle's viscosity r\ and rigidity fi . In the present geological 
epoch, the Maxwell time of the terrestrial mantle is about 500 years. For warmer mantles, 
it may be much shorter, taken the exponential temperature-dependence of the viscosity. 

Another characteristic time e ntering; the above expressions is the anelastic (Andrade) 



time t a . Referring the reader to lEfroimskyl (l2012al lbl) for details, we would mention that 



below some threshold frequency anelasticity seizes to play a major role in the internal friction, 
giving way to viscosity. Thus the mantle's behaviour becomes closer to that of a Maxwell 
body. Mathematically, this means that below the threshold the parameter r A increases 
rapidly as the frequency goes down. So only the first term in (130]) and the first term in (I3TT) 
survive, and we thus are left with the complex compliance of a Maxwell material. 



In our computations, we treated r A in the same way as in iMakarovl (120121 ) : we kept 
r 4 = r M at the frequencies above the threshold (which was set to be 1 yr -1 , just like in 
the solid Earth case). For frequencies lower than that, we set r A to grow exponentially 
with the decrease of the frequency, so the rheological model approached the Maxwell one in 
the low- frequency limit. Numerical simulation has demonstrated that the resulting capture 
probabilities are not very sensitive to how quickly the switch to the Maxwell model takes 
place. For more details, see Ibid. 

In a computer code, it is easier to divide both the numerator and denominator of (123]) 
by J 2 : 

_3 -A { X 

2(J - i) (k + a,) 2 + r 



e i(w ImM ) sin e/ (^ mw ) = — — Sgn(w ImM ) , (29) 



where 1Z and X are the dimensionless real and imaginary parts of the complex compliance 

~2 



K = 1 + ( X r A r a cos ( ^ ) r(a + 1) , (30) 
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1 



1-1 



sin 



a 7r 



(31) 



These dependencies were used also in iMakarovl (120121 ). to explore the spin-orbit dynamics of 
a Mercury-like planet. 



Appendix C. Capture probabilities 



T his sec tion offers a squeezed explanation of the capture theory developed by lGoldreich &: Peak 



( 119661 . 119681 ). 



As was mentioned in subsection 14.31 and explained at length in Appendix A above, a 
good approximation for the polar component of the tidal torque is furnished by expression 
([6]). The tidal-mode-dependent factors entering that expression can be exp ressed as functions 
of the spin rate 9 , see formula (jSJ) . Following iGoldreich fc Peald ( 119681 ) , in the vicinity of 
each particular resonance q' , i.e., for 9 being close to n (1 + q'/2) , it is instrumental to 
decompose the right-hand side of (jBJ) into two parts, one being the q — q' kink, another part 
being the bias. Constituted by the inputs from all the q ^ q' modes, the bias is smooth 
and virtually constant in the vicinity of the q ' resonance. 



The q = q' term is an odd function of w 2 2og' • While IGoldreich fc Peald ( 119681 ) assumed 
this term to be either a constant multiplied by Sgnu; 220 g ; (the constant-torque model) or 
a cons tant multiplied by the mode oo 220 q, itself (the linear in frequency model), Makarov 
(J2012J) endowed this term (in fact, all terms) with a realistic mode-dependence originating 
from the properties of actual rocks. 



It is convenient to introduce an auxiliary variable 



7 



29 - (2 + q')M 



(32) 



which vanishes when the long axis of the planet points toward the star at the perigee. Then 
the q = q' resonance will correspond to 7 = because 



1 



(2 + q')n + 29 



UJ, 



220( 



(33) 



Be mindful that our 7 c oincides with that e mploy ed by IMakarovl (120121 ) , and is twice the 
quantity 7 introduced in IGoldreich fc Peald (11968I ). 



In the absence of tidal friction, the equation of motion near the said resonance looks as 

M st ar 



C 7 + 3 (B - A) 



M. 



star 



M 



planet 



n G 2 oq'(e) sin 7 = 



(34) 



-31 - 



Multiplication thereof by 7 , with subsequent integration over time t , furnishes the first 
integral of motion, 



C — - 3 (B — A) Mstar n 2 G 20 q'(e) cos 7 = E' , (35) 

2 M st ar + Mpi ane t 

whose value depends on the initial conditions. The latter equation can be rewritten also as 

C ^ - 3 (B - A) Mstar n 2 G 20 cr(e) 2 cos 2 1 = E , (36) 

2 M star + M p i anet 2 



with E differing from E' by a 7— independent constant. As demonstrated in lGoldreich fc Peale 



( 119681 ). vanishing of the integral of motion E corresponds to the separatrix dividing rotations 



from librations. So the equation of this separatrix is 

3 (B - A) M t 1 1/2 



7 = 2 n 



L star 



G2oq'(e) 



C M star + M p i anet 
In terms of 7 , the q — q' term of the tidal torque reads as 



cos I . (37) 



W% = - KG z 220ql k 2 (j) sin|e 2 ( 7 )| Sgn (7) , (38) 

K being a positive constant. To write down the bias in the vicinity of this resonance, it is 
convenient to express an arbitrary tidal mode via the resonant one: 

u 220q = {2 + q)n - 29 = (2 + q')n - 26 + {q - q')n 

= - 27 + (q - q')n w (q - q')n , (39) 

where we recalled that 7 vanishes in the point of resonance. Taking ( 1391) into account, we 
can write the bias as 

V = K ^ G 2 20q k 2 (uj 220q ) sm\e 2 (u 220q )\ Sgn (u 220q ) 

= K ^2 G 22oqh((q - q')n) sin|e 2 ((g - q')n)\ Sgn (q - q') , (40) 

q+q' 



15 Our equations ([55)1 . and ([57)1 differ from their counterparts in iGoldreich fc Peald (|1968l ). The 

difference in numerical coefficients originates from our 7 being twice that used in Ibid. In our equations, we 



also keep the mass factor omitted in Ibid. 



16 



Recall that 



variable 7 



Mw 220(? ) sin|e 2 (w 220(? )| Sgn(w 220? ) 



is an odd function, wherefore the switch to the 



J 220q' 



generates a "minus" sign in 
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For a slowing-down planet, an estimate for the capture probability is derived from the 
consideration of two librations around the point of resonance 7 = — uj^ ng' = , i.e., of the 



l ast lib ration with a positive 7 and the first libration with a negative 7 . iGoldreich fc Peale 



( 119681 ) assumed that the energy offset from zero at the beginning of the last libration above 
the resonance is uniformly distributed between and AE = J (T)jdt . This assumption 
rendered them the following estimate for the probability of the capture: 

SE 

-T capt — ' V*1J 

with 8E being the total change of the kinetic energy at the end of the libration below the 
resonance. 

Thus, (T) 7 should be integrated over one cycle of libration, to obtain AE , and should 
be integrated over two librations symmetric around the resonance uJ 220 q, = , to obtain SE . 
As a result, the odd part of the tidal torque at q = q' doubles in the integration for SE , 
whereas the bias vanishes. Both these components are involved in the computation of AE . 

This makes capture probability look like 

p = 2 Z^wirddn 2 j: n wirt) d 1 = 2 

capt i: w [w(j) + v\d 1 ~ + 2 *v 1 + 2 tt vip w{i)6n ' 1 ] 



The in tegral in this equation can be evaluated if we further assume, following IGoldreich &: Peale 



( 119681 ). that in the vicinity of resonance the trajectory follows the singular separatrix solution 



BTj) which corresponds to vanishing of the integral of motion E . 
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Table 2. Default parameters of GJ 581 and planet d. 



Parameter Value 



R l-7-REarth 

Mpi ane t 7.1M Earth 

M star O.31M 

a 0.218 AU 

e 0.27 

(B-A)/C 5-10" 5 

Porb 67 d 

tm 50 yr 

H 0.8 • 10 11 kg m^ 1 s" 2 

a 0.2 
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Table 3. Symbol key 



Notation 



Description 



£ moment of inertia coefficient of the planet 

R radius of the planet 

T <TRI) triaxiality-caused torque acting on the planet 

T (TIDE) tidal torque acting on the planet 

M p i anet mass of the planet 

M st ar mass of the star 

a semimajor axis of the planet 

r ... instantaneous distance between the planet and the star 

/ true anomaly of the planet 

e orbital eccentricity 

M. mean anomaly of the planet 

C the maximal moment of inertia of the planet 

B the middle moment of inertia of the planet 

A the minimal moment of inertia of the planet 

n mean motion, i.e. ^ {M p i anet + M star )/a 3 

G gravitational constant, = 66468 m 3 kg -1 yr~ 2 

r M viscoelastic characteristic time (Maxwell time) 

t a anelastic characteristic time 

fj, unrelaxed rigidity modulus 

J unrelaxed compliance 

a the Andrade parameter 



